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Abstract 

I present the two-loop QCD corrections to the production of a massive quark- 
anti-quark pair in the massless quark-anti-quark annihilation channel. The result is 
obtained as a combination of a deep expansion in the mass around the high energy 
limit and of a numerical integration of a system of differential equations. The primary 
application of the outcome and developed methods is top quark pair production at 
the Large Hadron Collider. 

1 Introduction 

One of the most important goals of the Large Hadron Collider (LHC) is the measurement 
of top quark properties. This will be possible thanks to the abundant statistical samples 
reaching about 8 million pairs produced per year in the low luminosity phase. Besides the 
mass and decay parameters, the total cross section constitutes a primary observable. The 
experimental prospects of obtaining a measurement accuracy below ten percent for this 
quantity put a high demand on theoretical predictions of matching quality. At present, the 
known next-to-leading order corrections [1] have an error estimated from scale variation 
at about 12%. Soft gluon resummation [2-4], which has been an excellent tool for the 
Tevatron, and helped reduce the error to about 5%, is not safely applicable in the framework 
of the LHC. This is due on the one hand to the higher energy and on the other to the 
dominance of the gluon flux over the quark flux. Furthermore, the mentioned high statistics 
warrant the preparation of a Monte-Carlo generator of suitable precision, which cannot rely 
on resummed cross sections. 



In view of these facts it is necessary to provide a result for the next-to-next-to-leading 
order cross section, at best in a fully differential form. This requires four separate parts 
at the partonic level: the two-loop virtual corrections, the one-loop squared virtual cor- 
rections, the one-loop real-virtual corrections with an additional parton in the final state, 
and the tree-level real corrections with two additional partons in the final state. Out of 
these, the second part is known from [5,6], the third from [7], where the next-to- leading 
order corrections to the tt-|-jet corrections have been derived. Unfortunately, as part of a 
cross section calculation for top quark pair production this result still needs subtraction 
terms in order to allow for integration over the full phase space. Finally, while there is no 
result for the real radiation, the two-loop correction^ have been recently evaluated in the 
limit of small top quark mass [14, 15]. This result is applicable for highly energetic tops, 
for example when high pt cuts would be applied. The bulk of events comes, however, from 
the region much nearer to the partonic threshold. 

In this paper, I present a complete result for the two-loop corrections in the quark 
annihilation channel valid in the whole kinematically allowed region. It has to be stressed 
that obtaining an amplitude expressed in analytic form through some special functions 
seems out of reach in the nearest future. Since the LHC will soon become operational it is 
necessary to resort to semi- analytic /semi-numeric methods. The method adopted here is a 
combination of a deep expansion in the mass around the high energy limit, which contains 
the power corrections to the result of [14], and of numerical integration of differential 
equations. 

In the next section, I will first give a few definitions and then describe the power 
corrections. A detailed study of the numerical methods and the full result will follow in 
the last section of the main text. 



2 Power corrections 

The notation of this paper follows closely that of [14] . However, I reproduce all the neces- 
sary definitions for the convenience of the reader. 

The process under scrutiny is massless quark-anti-quark annihilation into a massive 
quark- ant i- quark pair 

ciiPi) + QiP2) QiP3,rn) + Q{p^,m) . (1) 
The amplitude can be described with the help of Mandelstam variables 

s= {pi+ P2 f, t= {pi- p^Y - u = {pi - pif - w?. (2) 

Notice that the mass subtraction in the definition of the t and u variables was irrelevant 
for the results of [14], because only logarithmic terms in the mass have been retained there. 
The advantage of this definition lies in the symmetric range of variation 
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"'^The case of massless quark production has been studied in a number of papers [8-13] 
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where (3 is the velocity 



The results will be parameterized by two dimensionless ratios 

rus = — , X = — . (5) 

s s 

The additional scale introduced by dimensional regularization, fi, has been kept in the 
results as unspecified, but for the plots and numerical values reproduced in the following 
jj, = m. The renormalization has been performed in the MS scheme with ni massless 
and Uh massive active flavors. The necessary constants are known in the literature with 
sufficient precision: the strong coupling renormalization to four- loop accuracy [16,17], and 
the mass and field renormalization of the heavy quark in the on-shell scheme to three-loop 
accuracy [18-20]. The renormalization of the light quark field, which is non- vanishing 
because of the presence of closed heavy quark loops has been explicitely given in [14] with 
two-loop accuracy. 

After expanding the amplitude in the strong coupling constant up to the second order 



\M) = 47ra^ 



(6) 



the interesting, two-loop term, contracted with the Born amplitude can be decomposed 
into color factors 

= 2Re(A^(°)|7W(')) = 2{N^ - 1) (7) 
X (^N^A + B + -^C + Nr^Di + Nn^D^ + ^Ei + "^E^ + n^^Fi + r^n^Fih + rj^^F^) . 

The leading behavior of the amphtude in the limit m — > 0, has been derived in [14] using two 
different approaches. The first is based on factorization properties of QCD amphtudes [21], 
and exploits a relation between the massless and massive cases. Unfortunately, it does not 
give a handle on mass corrections or the full mass dependence. The second approach 
is based on a direct evaluation of occurring integrals and is an evolution of a strategy 
developed for Bhabha scattering [22,23]. The procedure starts with a reduction of the 
integrals to masters with the help of the Laporta algorithm [24] , and subsequent expansion 
of the masters in the mass by passing through Mellin-Barnes representations [25,26]. The 
bulk of the work is performed by Mathematica packages MBrepresentation [27] and MB [28] 
together with further associated software. 

The derivation of the asymptotic behavior of MeUin-Barnes representations is performed 
recursively by closing the contours of integration and taking residues. It is obvious that 
arbitrarily high orders of expansion in the mass can be obtained. However, the coefficients 
will still be integrals requiring evaluation by summation of infinite series or by some other 
method. Previously, this last step has been completed with a combination of XSummer [29] 
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and the PSLQ algorithm [30]. It has to be stressed that every next order in e contains 
more integrals with a more complicated integrand structure. Fortunately, there exists an 
alternative approach based on differential equations [31,32]. 

Clearly, applying a derivative with respect to any invariant or the mass introduces 
higher powers of denominators and/or numerator structures. Furthermore, any set of 
integrals differing only by powers of denominators and numerators can be reduced to a 
smaller set of masters. In consequence one can write the following differential equation 
systems for the coefficients of the Laurent expansion of the master integrals 

d ^ , 



dm. 
d 
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dx^'^^"'^^ = ^J^jims,x)Ij{ms,x). (9) 

The Jacobian matrices, J^^ and J^, have rational function elements and it is implied that 
any master integral is a combination of the Ii{ms,x) functions 

Mi{ms,x,e) = ^e^/i^.(m^,x). (10) 

j=k 

The lowest power of e in the sum is fixed by the singularities of the integral and cannot 
exceed —4 at this level of perturbation theory, whereas the highest power is defined by the 
coefficient in the amplitude (there are spurious poles). In practice, the deepest expansion 
in e that occurred was down to order e'^ due to the particular choice of master integrals. 

The differential equations Eq. [8] allow, in principle, to fix the complete functional de- 
pendence of the master integrals. Unfortunately, the functions present in the solution are 
not known at present, and therefore a direct integration of the system has to be postponed. 
Nevertheless, the differential equations in can be solved by means of a series expansion, 
with boundaries given by the small mass limit as needed for the results of [14]. Following 
this idea, it is possible to derive arbitrarily deep expansions of the amplitude in the mass, 
and thus provide the power corrections, which were out of reach of the factorization ap- 
proach. Of course, the size of the intermediate expressions, combined with the available 
computing resources puts a natural cut-off on the highest power that can be computed. In 
fact, I have computed eleven terms of the series up to mj°. 

The results for the finite parts (in the e expansion) of the three purely bosonic con- 
tributions A, B and C, which are also the most involved as far as the computation is 
concerned are shown in Fig. [TJ The plots correspond to 90 degree scattering, i.e. x = 1/2. 
It is striking that the series do not obviously diverge, which is usually the case with this 
type of expansions. In fact, the series expansion for the leading color term, A, is at worst 
asymptotic at threshold, and can still be used to obtain an estimate accurate to a few 
percent at this point. On the other hand, the growth of the subleading color coefficients is 
indicative of the true behavior, but incorrect. As I will show in the next section, there is a 
true divergence due to the Coulomb singularity, which cannot be reproduced with a small 
mass expansion. 



4 




-100 



-150 \ < ~ < -=4 

2.5 5 7.5 10 

T\ = s/(4m) - 1 

Figure 1: Finite parts of the purely bosonic contributions to the two- loop amplitude as 
expansions around the small mass hmit at x — 1/2. The sohd hne represents the eleven 
terms of the series as derived for the present publication, the long dashed - ten terms of the 
series, the short dashed - six terms of the series, the dash-dotted - the leading behavior. 

The small mass expansion is, in reality, not an expansion in m^/s, but rather in 
max(m^/s, — m^/t, —w?/u). For small m and at the kinematical boundary corresponding 
to forward scattering 




A similar relation holds for —rn?/u for backward scattering. In consequence, the series will 
be asymptotic at best at the kinematic boundaries. 

In order to study the region of convergence and usability of the series, it is necessary to 
specify some criteria that could be applied without reference to any external result. In fact, 
if the amplitude were approximated with one permille accuracy, it would be sufficient for 
any foreseeable applications. Customarily, the error of an expansion is estimated by the size 
of the last term. For geometric series, this estimate is only correct (not underestimated) 
if the ratio of two subsequent terms does not exceed 1/2. In the latter case, eleven terms 



5 



0.75 




0.25 



T] = s/(4m ) - 1 



Figure 2: Convergence regions of the small mass expansion of the two-loop amplitude. 
The grey area represents the region where the series is accurate to one permille, whereas 
the black area the region where the leading term of the series is accurate to one percent. 
The dashed lines delimit the kinematically allowed region, whereas the short dashed lines 
inside the grey area would correspond to the convergence region derived according to Eq. 
without the last condition. 



of the series (as in the case of the present result for the amplitude) provide indeed an 
approximation exact to one permille. In the case of amplitudes, the series is obviously not 
geometric, but conditions inspired by these arguments can be imposed. Let the amplitude 
be written as 



^(0,2) 



10 



i=0 



(12) 



A relatively conservative heuristic test for the one permille convergence of the result is 
given by the following conditions 



aiQm 



10 



< 10"^ /\ 
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Og 




aiQ m 



10 



< 10' 



(13) 



which are applied at a given (m<(, x) point during a scan starting from the median x = 1/2, 
which has always the best convergence for a given value of m^. The last condition in 
Eq. [13] deserves further explanation. It turns out that for relatively small values of m^, 
the logarithmic terms in lead to slight violations of the remaining relations, but the 
last term of the series is still tiny. In this case, it is highly improbable that the sum of 
the missing terms would amount to more than a permille correction. Without the last 
test the region of permille convergence is, therefore, unrealistically restricted. The region 
resulting from the application of Eq. [13] is shown in Fig. [2l The visible discontinuities of 
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the boundary are due precisely to the nature of the criterion (and to a lesser extent to the 
discretized scan). This figure shows also the region, where the leading term of the series 
agrees with the full result to one percent. No special criteria are needed here, of course, 
since eleven terms of the series expansion provide a sufficient approximation to the exact 
result for the purpose of determining this region. 

3 Numerical solution 

Since the series expansion does not satisfactorily approximate the result over the whole 
range of variation of kinematic parameters, it is only natural to try to solve the differ- 
ential equations for the master integrals numerically and thus obtain a purely numerical 
description of the amplitude. This idea has originally been put forward in [33] for mas- 
ter integrals corresponding to the two-loop sunrise graph without, however, relating to a 
concrete physical problem. In fact, to the best of my knowledge, the only applications to 
physical processes have been attempted in [34,35]. The problem at hand pushes the diffi- 
culty level substantially further, and requires, therefore, a careful assessment of feasibility. 
Before 1 show that high precision numerical results may indeed be obtained at all relevant 
points of phase space, there are several issues that have to be clarified, as it can hardly be 
overstressed that the right choice of numerical algorithms will make the difference between 
success and failure. 

3.1 Implementation 

At first, it is necessary to determine the type of differential equation system under consider- 
ation. It turns out that different methods are available for stiff and non-stiff problems. The 
distinction between the two is somewhat fuzzy in the professional literature, but there is 
agreement that stiff problems involve exponentially decaying solutions. The criterion used 
in practice is the existence of large negative eigenvalues of the Jacobian. In the present case 
it is, however, easier to use heuristic arguments instead of performing a numerical analysis. 
In fact, experience accumulated in numerous higher order calculations shows that expo- 
nentially decaying components would be rather unusual. I will, therefore, assume without 
further consideration that the system is non- stiff. In this case, there is still a large number 
of algorithms available. However, because the solutions, i.e. the master integrals, must 
be very smooth (we remain above all thresholds) and high precision will be requested, a 
variable coefficient multistep method [36] is expected to be most efficient [37]. Indeed, 
this kind of methods is based on polynomial interpolation/extrapolation with polynomials 
of order up to 12, which is a guarantee of very fast convergence under the assumption of 
suitable smoothness (if higher order derivatives were large, the errors would obviously grow 
severely with the order). 

The next choice concerns the treatment of master integrals, as they can be considered 
to be either complex or two-component real functions (after decomposition into real and 
imaginary parts). Clearly, working with complex functions reduces the size of the Jacobian, 
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which may give substantial improvements in the running time. Furthermore, the real 
function approach suffers from the fact that the imaginary parts of many integrals vanish 
for real arguments, which poses problems as far as error control is concerned. Indeed, in 
such cases it is impossible to use a relative error criterion and absolute errors must be used. 
The size of the latter can only be determined in conjunction with the final amplitude in 
order to have control over the precision of the outcome. There seems to be only one, albeit 
very strong, argument in favor of real functions, namely that most of the available software 
works with real numbers exclusively. A quick glance at the literature of the subject shows, 
that writing a code from scratch should be avoided. Fortunately, one of the most advanced 
software packages implementing the variable coefficient multistep method [38] has recently 
been translated to complex arguments. 

Closely connected to the choice made above is the problem of error control. Customarily, 
working with complex functions implies that the error is given by the modulus of the 
difference between the exact value and the approximation. In the present case, however, 
we are only interested in the real part of the amplitude, hence the imaginary components 
of the master integrals will be discarded. In consequence, we need to control the error 
of the real part and not that of the modulus. Fortunately, unless the imaginary part is 
much larger than the real, the two are not much different. It turns out that in the present 
calculation, only about 6% of the evaluated phase space points involved an integral, for 
which the imaginary part was more than lO'^ larger than the real part. Therefore, for 
simplicity reasons, the traditional error estimate has been used in the following. Notice 
also that the imaginary parts could have been discarded from the start, since the system 
of differential equations is linear and we do not cross any singularities. 

After settling the implementation questions, it is necessary to decide on the position 
of the boundary conditions. In the original publication [33], it has been proposed to start 
from a threshold or a pseudo-threshold, since the values of the integrals at these points were 
known. This approach has the drawback that these points are at the same time singularities 
of the differential equations, which requires slight modifications of the algorithm and leads 
invariably to a substantial loss of precision when evolving further from the boundary. Here, 
I use the series expansion of the previous section to compute the values of the integrals to 
very high precision. In fact at 

= 5 X 10"^ ^ ^ i' ^^^^ 

the relative error estimated by the size of the last term (very conservative) does never 
exceed 10^^^. The second condition in Eq. [T3] deserves explanation, because the median 
point X = 1/2 would have led to better convergence. Unfortunately, we will see later that 
it is also a singular point of the system and thus cannot be used. The choice taken results 
in the loss of about two digits and is compensated by a twice smaller value of m^. 

Let me now turn to the discussion of the contour of integration. It is clear [33] that 
evolving along the real axis should be avoided, in order not to stumble on the singularities. 
Since we are working with complex functions anyway, it is not a problem to deform the 
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Figure 3: Evolution contours for the numerical integration of the system of differential 
equations in and x. 

contour into the complex plane. In general, this deformation is only restricted by causalitjl^. 
In the present case, however, we always remain above thresholds, which means that the 
Riemann sheet has already been chosen when fixing the boundary conditions, and any 
curve will be appropriate. To take full advantage of the multistep algorithm for integration, 
which depends on several previous values of the system, it would be desirable to use a single 
smooth curve to reach a given point starting from the boundary. Unfortunately, we need 
an integration in two variables, and experimentation has shown that it is more efficient to 
perform two separate evolutions. For each of these, I will use an elliptic contour, with a 
user specified eccentricity as depicted in Fig. [31 The latter freedom allows for a relatively 
easy estimate of the final global error, by computing the desired amplitude with several 
different contours. It is also interesting to note that in practice the convergence turns out 
to be faster for more circular contours. 

A rather unpleasant feature of the system of differential equations at hand is the pres- 
ence of a number of singularities in the Jacobians. They are summarized in Tab. [H It is 
crucial that aside from thresholds the solution is regular at these points. Therefore, in the 
case of the four different singularities, which occur inside the kinematically allowed region 
it is sufficient to resort to interpolation. The problem is more acute for the singularities at 
the boundary (forward/backward scattering). These approach the true branching points 
at the thresholds when m — 0, and therefore interpolation is only efficient for moderate 
values of m, when the necessary points outside the kinematically allowed region are not 
trapped too close to the branching points. Otherwise, extrapolation is necessary. Notice 
that the concept of "dangerous distance to a singularity" requires specification. In fact, it 
is defined by the available numerical precision and the strength of the singularity (power of 
the singular polynomial in the denominator of a coefficient). Finally, one has to remember 
that the solution for the master integrals will be input into the expression for the am- 

^for example, if the evolution were performed in the Mandelstam s variable, the contour should be in 
the upper half-plane, when approaching the cut. 
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rris = l/2x 




yes 




rus = 1/2 (1 - x) 




yes 




rris = 1/2 (1 - x^) 








ms = -1/2(1 -x)2 









Table 1: Complete list of singularities of the Jacobians, J and J , of the system of 
differential equations. The singularities occur in both systems of differential equations 
in rris and x apart from the point rris = "1/4? which is present only in the differential 
equations in m^. The table indicates in addition the presence of a branching point at 
a given singularity (a blank entry denotes a regular point of the solution), and specifies 
whether a singularity occurs within the kinematically allowed region (a blank entry denotes 
a point outside the allowed region) . 

plitude, which may lead to further cancellations, and hence instabilities. For the present 
calculation, I have adopted extended precision (quadruple) in the numerics to overcome 
this problem. 

3.2 Efficiency and numerical stability 

In order to illustrate the efficiency of the approach, I show in Tab. [2] different timings and 
other related informations for the complete solution at the point rris = -2 and x = .45. 
Since the evolution is performed separately first in rris and then in x, I require two more 
digits of local precision in the first step, so that the estimate of the error will be dominated 
by the second. Note that in the first evolution the x value is fixed from the very beginning 
as specified in the boundary condition Eq. [TH This value has been hardcoded in the 
Jacobian, which not only leads to a much faster evaluation, but also to less severe numerical 
instabilities in the coefficients themselves. The latter feature is actually crucial in reaching 
the higher precision in the first evolution. Returning to the error estimate, it has to be 
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leading eolor 


full color 


number of masters 








36 






145 




number of functions 








155 






595 




precision 


quadruple 


double 


quadruple 


double 








evolution in m 










requested local error 






10-20 


10-12 


10-12 


10-20 


10-12 


10-12 


contour deformation Snis 


0.1 


0.1 


0.1 


0.1 


0.1 


0.1 


number of steps taken 




2319 


618 


534 


2932 


777 


1302 


Jacobian evaluation time [ms] 


3.4 


3.4 


0.2 


37 


37 


4.9 


evolution in x 


requested local error 








10-10 


10-10 


10-18 


10-10 


10-10 


contour deformation 5x 




0.1 


0.1 


0.1 


0.1 


0.1 


0.1 


number of steps taken 




545 


139 


139 


739 


174 


432 


Jacobian evaluation time [ms] 


8.3 


8.3 


0.4 


150 


150 


17 


total evaluation time 


s 




49 


13 


< 1 


957 


243 


26 



Table 2: Timing and efficiency information for the numerical integration of the system of 
differential equations to the point — .2, x — .45. The numbers have been obtained 
on a 2GHz Intel Core 2 Duo system, after compilation with the Intel Fortran compiler. 
Quadruple precision is an option of the compiler. 



stressed that in numerical integration of systems of differential equations, the precision is 
specified locally, i.e. one requires that the error of the approximation does not exceed a 
certain value at every step. Therefore, the global relative error is estimated by the product 
of the number of steps taken and the requested local error. In practice, the number of 
steps in the evolution never exceeded ten thousand, and therefore the final precision 
for a local relative error of IO-20 (in quadruple precision) was roughly sixteen digits. The 
number of steps needed in the second evolution is usually smaller, because the starting 
point is far from any singularities. In consequence for the requested local error of 10"^^, 
the final global error should not exceed about lO-i^. There is one additional source of 
error connected to roundoff and numerical cancellations in the coefficients. Its presence 
is visible in the table, when comparing the solution in double precision and quadruple 
precision for the same requested local relative error. In the case of the x evolution of the 
full system, the number of steps is larger in double precision, precisely because the error 
estimates are not satisfied due to random roundoff errors. The software tries to reduce the 
step size until the error estimate satisfies the bound, which eventually happens because 
the random variations around the true value must, sooner or later, turn close to it. Let 
me finally comment on the running time. Clearly, if only the leading color coefficient were 
needed with moderate precision, a value at a single phase space point could be obtained 
within less than a second. For the full color structure in quadruple precision, as much as 
fifteen minutes are needed. Although this does not allow for direct implementation in a 
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0.22625 


-0.06808683395 


-18.00716652 


6.302454931 


3.524044913 


Di 




-0.22625 


0.2605057339 


-0.7250180282 


-1.935417247 


Dh 






0.5623350684 


0.1045606449 


-1.704747998 


El 




0.22625 


-0.3323207300 


7.904121951 


2.848697837 


Eh 






-0.5623350684 


4.528240788 


12.73232424 


Fi 










-1.984228442 


Fih 










-2.442562819 


Eh 










-0.07924540546 



Table 3: Values of the color coefficients of the two- loop amplitude at the point = 
.2, X = .45 rounded at 10 digits precision (the given digits are unaffected by numerical 
uncertainties). 



Monte-Carlo generator, the functions are smooth enough to be interpolated starting from 
a grid of values. The efficiency is by far sufficient to obtain dense grids. Therefore, a grid 
of numerical values for moderate rris together with the series expansion of the previous 
section for small rris is a complete solution to the problem of evaluation of the two-loop 
amplitudes for the production of a heavy quark-anti-quark pair in massless quark-anti- 
quark annihilation. In Tab. [3l I give the values of all the color coefficients with ten digits 
precision at the point = .2, x = .45, which is well outside the region of convergence 
of the series expansion. The actual precision of the result estimated by the variation with 
respect to the change of the integration contour in x was roughly fourteen digits. 

In the course of preparation of the present work, I derived the necessary grids mentioned 
above. For simplicity the singularities at the kinematic boundaries (forward/backward 
scattering) have been avoided by keeping a distance of 10~^ in x, which is unnoticeable on 
the plots, but can be improved upon in the future. In this respect, the actual needs will 
only be apparent once the virtual corrections will be combined with real radiation. The 
range of variation in rus was chosen such that the distance to the threshold parameterized 
by 

4m^ 

was at least 10"^. This is safely sufficient for any practical applications. The plots of the 
finite parts of the purely bosonic corrections can be found in Fig. |H whereas those for the 
contributions containing a single closed fermionic loop in Fig. [51 The interesting feature is 
the large variation of the subleading color contributions, when nearing the threshold. Of 
course, this is due to the Coulomb singularity. In particular, the Cq coefficient behaves 
like l/jS^, which cannot be compensated by phase space integration. This leads to a true 
divergence, to be taken care of by resummation in the complete analysis. 
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Figure 4: Finite parts of the bosonic contributions to the two- loop amphtude. 

Finally, let me point two possible improvements of the implementation. The first one is 
related to the numerical precision. Clearly, different choices of the basis functions may be 
used to milden the strength of the singularities. For example, for two functions, which suffer 
from large cancellations, one could introduce a mixture such that one of the functions is 
small, while the other is large, but contributes with a small coefficient. The second possible 
improvement is at a lower level and concerns the time of evaluation. Since the computation 
is done in quadruple precision, one could try different libraries instead of the compiler's 
built-in routines. In fact, the implementation of [39] has proven to be about three times 
faster on some problems. Of course, this is of lesser importance than precision, since even 
times of the order of five minutes per point would still not allow to perform the integration 
in real time within the framework of a Monte Carlo generator. 

4 Conclusions 

In this paper, 1 have presented a complete solution to the problem of evaluation of the 
two-loop amplitude for heavy quark production in light quark annihilation. The result. 
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Figure 5: Finite parts of the single fermionic contributions to the two- loop amplitude. 

a numeric approximation obtained by a combination of a small mass expansion and in- 
tegration of a system of differential equations, is satisfactory from the point of view of 
applications. Aiming at a complete description of the top quark pair production cross 
section at the next-to-next-to-leading order at the LHC, the next steps to perform will be 
the evaluation of the two-loop amplitude in gluon fusion and the computation of the real 
radiation contributions. The mixed real-virtual contributions are in principle known, but 
have to be supplied with suitable subtraction terms in order to allow for integration over 
the full phase space. 

The result for the series expansion of the amplitude, as well as the grid of values obtained 
by numerical integration of differential equations (with all numbers rounded at 5 digits) 
are available in the form of Mathematica files, qqQQ2Lseries .m and qqQQ2Lnumeric .m 
respectively, attached to the source of the paper on the preprint server |http : / /arXiv . org| 
but can also be obtained from the author upon request. 
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